arXiv:1503.06226vl [physics.flu-dyn] 20 Mar 2015 


Laminar-turbulent patterning in wall-bounded 
shear flows: a Galerkin model 

K. Seshasayanan 1 and P. Manneville 2 

1 Laboratoire de Physique Statistique, CNRS UMR 8550, 

/ 

Ecole Normale Superieure, F-75005 Paris, France 

2 Laboratoire d’Hydrodynamique, CNRS UMR7646, 

/ 

Ecole Polytechnique, F-91128, Palaiseau, France 
to appear in Fluid Dyn. Res. 


Abstract 

On its way to turbulence, plane Couette flow - the flow between counter-translating 
parallel plates - displays a puzzling steady oblique laminar-turbulent pattern. We ap¬ 
proach this problem via Galerkin modelling of the Navier-Stokes equations. The wall- 
normal dependence of the hydrodynamic field is treated by means of expansions on func¬ 
tional bases fitting the boundary conditions exactly. This yields a set of partial differential 
equations for the spatiotemporal dynamics in the plane of the flow. Truncating this set 
beyond lowest nontrivial order is numerically shown to produce the expected pattern, 
therefore improving over what was obtained at cruder effective wall-normal resolution. 
Perspectives opened by the approach are discussed. 
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1 Context 

Turbulent flows display transport properties strongly enhanced with respect to those of laminar 
flows, a feature that has particularly important consequences in configurations of engineering 
interest. Understanding how a given laminar flow becomes turbulent or a turbulent flow decays 
to laminar is therefore of great interest, both conceptual and practical. In this respect, the 
case of wall-bonnded flows is of utmost concern since the transition can be direct, without 
the intermediate steps observed, e.g. in free shear flows (Huerre & Rossi 1998, Manneville 
2015). This direct transition is a result of the local stability of laminar flow in competition 
with nontrivial solutions to the Navier-Stokes equation (NSE for short) in the range of control 
parameter where the transition effectively takes place. As analysed by Waleffe (1997), the 
mechanism by which such nontrivial solutions exist, the self-sustainment process (SSP), is now 
thought to be well understood, but the laminar-turbulent coexistence still raises important 
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questions. In Hagen-Poiseuille flow (HPF), the flow under pressure gradient through a circular 
pipe, the transition takes place when turbulent puffs , the nontrivial states alluded to above, 
split and propagate turbulence before they have time to decay, a scenario well-reproduced 
by a reaction-diffusion model introduced by Barkley (2011a). In its own transitional range, 
plane Couette flow (PCF), the simple shear flow developing between counter-translating plates, 
experiences laminar-turbulent coexistence in the form of steady oblique bands (Prigent et al. 
2002, Duget et al. 2010, Barkley & Tuckerman 2005). The Reynolds number, the relevant 
control parameter, is here defined as R = Uh/v where U is the speed of the plates, h the half 
gap width between them, and v the kinematic viscosity of the fluid. The bands are observed 
for R g < R < R t . Below the global stability threshold R g turbulence is only transient, in the 
form of finite-lifetime spots, and the laminar base flow is always recovered after the spots have 
decayed. Beyond the upper threshold R t turbulence is essentially featureless , i.e. uniform. A 
model, also of reaction-diffusion type, was proposed by one of us (Manneville 2012) to account 
for this pattern formation, in which a Turing mechanism was proposed to be responsible for 
the bands when R is decreased below R t . Such explicative models are analogical in essence. 
Trying to support them directly from a reliable simplification of the primitive problem in order 
find out the physical mechanisms behind laminar-turbulent coexistence is the actual purpose 
of the work presented here. 

In the transitional range, the nontrivial solutions appear to be strongly coherent at the scale 
of the distance to the wall, pipe diameter (Hof et al. 2004) or gap between plates (Bottin et 
al. 1998). This feature can be incorporated in the sought-for models using Galerkin methods 
that project the solutions and their dynamics on well-chosen functional bases (Finlayson 1972). 
Such a method was used by Waleffe to build a dynamical system implementing the SSP for 
PCF with stress-free boundary conditions directly from the NSE (Waleffe 1997). His model 
is a system of ordinary differential equations governing the amplitude of velocity components 
involved in the SSP upon assuming full coherence at the scale of a Minimal Flow Unit (Jimenez 
& Moin 1991), with size of the order of the distance between the walls. Though valuable to 
discuss the SSP, this assumption is not appropriate to study the extended systems of interest, 
long pipes or wide channels. 

The method can however be adapted to such cases for which the variables have to be held 
amplitudes governed by partial differential equations still involving spatial coordinates rather 
than scalars functions of time satisfying ordinary differential systems. Basically, coherence in 
the flow is taken into account by projecting the hydrodynamic variables onto a limited number 
of wall-normal modes with corresponding amplitudes depending on the remaining coordinates, 
axial for pipe flow, in-plane for PCF. In spirit, this modelling approach can be considered as the 
analytical implementation of a direct numerical simulation (DNS) scheme in which the wall- 
normal spectral resolution would be varied in a controlled fashion. When practiced in a strictly 
numerical context (Manneville & Rolland 2011), this strategy showed that laminar-turbulent 
bands in transitional PCF are remarkably robust since they are preserved upon drastically 
reducing the number of Chebyshev polynomial used to represent the wall-normal dependence 
of the flow. As the resolution was decreased, the quantitative price to be paid was a progressive 
narrowing of the Reynolds number range where the bands were observed, accompanied by a 
downwards shift of that range explained by an aborted energy transfer towards small scales, 
whereas the main properties of the pattern were preserved qualitatively speaking. 

Performing more work “by hand” would yield equations that would be more cumbersome 
than the NSE but would encode wall-normal coherence at the moderate Reynolds numbers 
of interest in a crucial way. Eliminating supposedly less relevant degrees of freedom, we can 
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hope for better physical understanding and for higher numerical efficiency since the physical 
dimension of the problem would then be reduced from three to two. The limit would of course 
be that the obtained equations be still manageable. Here, we present an extension of a previous 
Galerkin model by Lagha and one of us (Lagha & Manneville 2007a) that pursues this agenda. 
In that study, the model was truncated at lowest significant order and retained only three fields. 
It displayed most of the expected qualitative properties except for the presence of any organised 
laminar-turbulent coexistence in wide domains (Manneville 2009). Taking into consideration 
the previously mentioned simulation results at reduced resolution (Manneville & Rolland 2011) 
we surmised that this deficiency would be corrected by truncating the expansion at a higher 
level. With four more fields, numerical simulations display the expected patterns as described 
in section [3j therefore validating the approach sketched in section [2} Perspectives opened by 
this approach will be discussed in section [4j The explicit expression of the model is given in [Aj 

2 Model 

The derivation follows previous work in (Lagha & Manneville 2007a) with the difference that, in 
order to avoid difficulties in the treatment of the pressure field, the NSE governing the departure 
from laminar flow is now written in a velocity-vorticity formulation as described in (Schmid & 
Henningson 2001), p.l55ff, i.e. the (nonlinear) Orr-Sommerfeld equation for the wall-normal 
velocity component v : 

0 d t + u h d x )V 2 v - ul{p x v + Af v = vV A v , (1) 

and the Squire equation for the wall-normal vorticity component (, = d z u — d x w, where u and 
w are the streamwise ( x ) and spanwise (^) velocity components, respectively: 

(d t + u h d x )( + u' h d z v + Af c = uV 2 C. (2) 

In these general equations the base flow is Vb = u^(y)e x . When dealing with PCF, using the 
half-gap width h as length unit and h/U as time unit, with U the speed of the plates driving 
the flow at y — ±1 we have = y. In that system of units the Reynolds number R is 

just numerically equal to l/i/, i.e. the inverse of the kinematic viscosity of the fluid. Primes 
denote the differentiation with respect to the wall-normal coordinate ( y ), hence u' h = 1 and 
u'^ = 0. The nonlinear terms AF v and Af( are complicated, formally quadratic, expressions of the 
velocity components and their derivatives that can be found in (Schmid & Henningson 2001). 
It will turn out convenient to use a poloidal-toroidal decomposition of the hydrodynamic fields 
by introducing a velocity potential 0 and a stream function 0 such that: 

v = -A 0, (3) 

u = -d z ' 0 + d xy 4 >, w = d x i4 + d zy 4> , C = ~M >, (4) 

A denoting the in-plane Laplacian d xx + d zz . 

The Galerkin approach used in (Lagha & Manneville 2007a) separates the in-plane space 
dependence of the hydrodynamic field from its wall-normal dependence by expanding it onto a 
polynomial basis in y, yielding amplitudes functions of x, z and time t. The no-slip boundary 
conditions to be fulfilled read ■u = r? = u’ = () = '0 = 0 = Oat 2 / = ±1 and, from the continuity 
condition d x u + d y v + d z w = 0, d y v = d y (j) = 0 at y — ±1. The basis functions are chosen so as 
to fulfil these boundary conditions exactly. The functions chosen for u , w, and 0 are in the 
form fi(y ) = (1 — y 2 )Ri(y), i = 0,... ,i max ; the Ri are polynomials of degree i, and i max is the 
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Basis functions for ’Q Basis functions for v 




Figure 1: Basis functions for the wall-normal vorticity ( (left) and the wall-normal velocity v. 
Functions used in our model are displayed with thick lines; / 0 : black, dotted; (/i,<?i): green, 
continuous; (/ 2 , 5 ' 2 ): red, dashed; (./s, . 9 , 3 ): blue, dash-dotted. Higher-order functions are shown 
with thin lines. The work in (Lagha & Manneville 2007a,b) made use of {fo,fi,gi} only. 


truncation order. For v and the functions are taken as gi(i ) = (1 — y 2 ) 2 Si(y), i = 1,..., 
the Si are polynomials of degree i — 1 for consistency with the continuity condition at given 
j m ax' The bases {/*} and {g*} are separately made orthonormal via a standard Gram-Schmidt 
procedure using the canonical scalar product (r|s) = \ r(y)s(y)dy. Basis functions are 
shown in figure [I] from which it is clearly understood how (z) the resolution close to the plates 
is improved by increasing the truncation order, and (ii) the profiles chosen for v incorporate 
the boundary condition d y v(x, ±1, z, t) = 0. The analytic expressions of basis functions up to 
j m ax = 5 are given in [AT As pointed out by Rolland (2012) in Appendix B of his PhD thesis, 
the chosen basis {fi,gi} is related to Jacobi polynomials of alternate possible use in standard 
spectral methods for the NSE (Canuto et al. 2007). 

According to the standard Galerkin procedure (Finlayson 1972), the expansions {v, 0} = 
J2 i {V i ,$ i }(x,z,t)g i (y) and = Yji{ U iiW il Z i ^ i }(x,z,t)f i (y) are inserted in the 

equations which are then projected onto the relevant bases, Eq. [ljfor v onto {g«}, and Eq. [ 2 ] for 
f onto {fi}- The concrete derivation is straightforward and can be automated once the order 
of truncation z max has been fixed. The formal expression of the model reads: 


{(IA + A) d t + (IA + d x — u (IA 2 + 2 AA + P)} A<J> = N (y) , (5) 

{(I d t + B<9 X ) - v (IA + P)} AT + B(9 2 A<F = N (z) . ( 6 ) 

In these expressions <f> and T respectively stand for arrays {$!,..., <Fj max )} t and 

{T 0 ,..., 'I f j rnax )} t , superscript ‘t’ denoting transposition. I is the identity matrix of order z max 
in Eq. [ 5 ] and z max + 1 in Eq. [6} I is a square but non-diagonal matrix of order z max + 1, playing 
a role similar to I. All other matrices, A, A, B, B, P, and P are either square or rectangu¬ 
lar, with coefficients straightforwardly obtained by integration over [—1,1] of the appropriate 
products of fi, gj and their derivatives. About one half of the possible combinations cancel 
due to parity considerations. Remarkably enough, matrices A, P and P are diagonally dom¬ 
inant [(i,i) (i,i + k)\ and the absolute value of the diagonal terms increases rapidly with 
the position of the coefficient <C (i + l,i + 1)], which suggests possible simplifications 

in the equations governing the dynamics of the field amplitudes. Finally, -> and are 
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complicated, formally quadratic expressions of the (C/j, Vi, W*)’s that have to be derived from 
the ($j, 'hj)’s introduced upon elimination of the pressure. Their explicit expressions are given 

in [Ol 

Equations ([5]-[6]) only involve A<3> and AT, which implies that some care is needed when 
dealing with spatially averaged terms corresponding to Fourier modes at (k x , k z ) = (0, 0). This 
is the price to pay for having used the velocity-vorticity formulation that avoid the explicit 
treatment of the pressure held (Schmid & Henningson 2001). It is then convenient to identify 
the uniform contributions to u and w explicitly by assuming: 

u = u — d z 0 + d xy 4 >, w = w + <9 X 0 + d zy 4 >, 

with u and w still function of y and t but independent of x and z, while 0 refers to the (x, z)- 
varying part of 0. Notations being unambiguous, the tilde will be dropped in the following. 
On general grounds, the mean how components u and w are governed by 

d t u — v(u)" = — (uv )', d t w — u(w)" — — (urn )', (7) 

where the overline means averaging over the in-plane coordinates. In the model, this is treated 
by expanding u and w onto basis {/)}. From the continuity equation we get: 

U = U - d z y + C d x $ , W = W + + C d z § (8) 

where U and W stand for arrays {U o, ■ ■ ■ U t max } T ' and {Wo,... } r ' while matrix C arises 

from the projection of d y v onto the basis used to expand u and w in the continuity equation. 
Upon projection, equations [7] read: 

ld t U - vWU = , I d t W - vFW = Nq Z) (9) 

where and No“' ) are the projections of the (x, z) spatially averaged nonlinear terms in 
equations [7J The model is now complete and ready for use. 


3 Validation 

By construction, the model possesses all the properties requested to account for the transitional 
regime of PCF: it can be checked that laminar how is linearly stable for all Reynolds numbers, 
despite transient energy growth linked to lift-up, and that its nonlinearities redistribute but 
conserve the kinetic energy contained in finite amplitude perturbations. A numerical solver 
was developed in order to examine whether bands can be recovered beyond lowest nontrivial 
truncation order i max = 1. Wanting to add higher modes of both parities, we chose i max = 3, 
i.e. 7 helds: Tj, i — (0 : 3), and $«, i — (1 : 3). In-plane space dependence was handled using 
a Fourier pseudo-spectral scheme that gets rid of aliasing via the usual 3/2 rule (Canuto et 
al. 2007), i.e. in directions {x,z}, the numbers of evolving modes are N{ XtZ \ and nonlinear 
terms are evaluated via back-and-fro FFTs with solutions reconstructed on |A/ x .. 2 } points. 
Time marching was treated by formally rewriting the initial problem d t X = CX + M{X) as 
<9i[exp(— tC)X] = exp(— t£)J\f(X) and solving the new problem using a Runge-Kutta scheme 
of order 4. 

In parallel to the study in (Manneville & Rolland 2011) dealing with DNS at reduced wall- 
normal resolution, a hrst numerical experiment was devoted to the recovery of the featureless 
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Figure 2: Fourier spectral power of the streamwise perturbation velocity component u in the 
mid-plane y = 0 as functions of wave-numbers n z for n x = 0 (left) and n x for n z = 0 (right). The 
corresponding Fourier wave-vectors read ky x>z y = 2nny X!Z y/Ly X}Z y, where Ly x>z y are the stream- 
wise (x) and spanwise (z) dimensions of the computational domain. The curves correspond to 
the different resolutions studied, in the form N x * N z , where Nr XjZ y are the maximum running 
wave-numbers. 

turbulent state belonging to the nontrivial branch at high R in a domain of size L x xL z = 32x32, 
aiming at an optimisation of the in-plane numerical resolution in view of the reliable simulation 
of domains as wide as possible at the cheapest possible computational cost. Results shown in 
figure [2] display the power spectra of the streamwise component of the perturbation velocity u as 
a function of wave-numbers n z for n x = 0 (left) and of n x for n z = 0 (right). Normalisation by 
the total number of modes N X N Z makes the curves corresponding to the different resolutions lie 
on top of each other. In the left panel, the peak generated by the spanwise statistical periodicity 
of streaks and streamwise vortices is clearly identified for all the resolutions considered but more 
pronounced for N x x N z = 128 x 128 than for 32 x 32. This corresponds to N XyZ /L XyZ = 4 and 
1, with effective space steps S x>z = 0.25 or 1, respectively, to be compared to the period of 
the streaks A~ ~ L z /n s t r with n s tr ~ 7, hence about A^ ~ 4.6 in reasonable agreement with 
known results. As seen in the right panel, the streamwise correlations decrease in a monotonic 
way as expected from the discussion in (Philip & Manneville 2011) where size effects on the 
temporal vs. spatiotemporal character of the dynamics was scrutinised. The subsequent study 
is restricted to the configurations mentioned in the table below, with the resolutions found 
acceptable from the previous experiment: 


Lx 

L z 

N x / L x 

N z /L z 

108 

48 

2 

4 

128 

84 

2 

4 

680 

340 

1 

1 


Our main result is that, in all cases, steady oblique patterns of alternately laminar and 
turbulent domains were observed in a limited range of Reynolds numbers, between R g zz 150 
and Rt ~ 159. Figure [3] (top) for R = 151 illustrates the two different possible orientations of 
a single band pattern in the domain 108 x 48. Orientation fluctuations are known to exist in 
DNSs at such an intermediate size. They are also present in the model as seen in the bottom 
panel showing the alternative dominance of modes (1,+1) and (1, —1), while the other modes 
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Figure 3: Patterning in a domain of size 108 x 48 at R = 151. Top: the two different orien¬ 
tations (1,+1) (left) and (1,-1) (right) with perturbation energy field averaged over the gap 
in grey levels, black = laminar, white = largest local energy. Bottom: Orientation fluctuations 
evidenced by the spectral power in modes with wave-numbers (n x , n z ), n x = 1, 2, n z = ±1, ±2. 

are less intense. Here fluctuations seem more important than in DNSs, with briefer episodes of 
well-formed pattern and a much smaller signal-to-noise ratio; compare with Fig. 3 of (Rolland 
& Manneville 2011). 

Larger domains can accommodate more bands as seen in figure [f] (left) for a 128 x 84 domain. 
The average of the turbulent energy over the volume V of the domain, E t — y f y |(u 2 + v 2 + 
w 2 ) dx d y dz, has been measured through the whole transitional range. The bifurcation diagram 
displayed in figure [4] (right) is again as expected, however the occurrence of large-scale laminar- 
turbulent coexistence in the form of oblique bands, easily detected visually and permitting the 
identification of R g ~ 150 and Rt ~ 159, leaves a weaker signature on the variation of E t with 
R than in DNSs for which a marked break at R t and a linear decrease below were observed. 
Here the smoother variation of E t (R) and the absence of clear-cut change at R t between the 
band regime (open circles) and uniform turbulence (filled circles) are presumably again a direct 
consequence of the higher level of fluctuations. Below R = 150, turbulence is only transient but 
a mean energy, roughly constant before the decay stage, can still be measured (open square). 
In principle R g should be located using a statistical study in line with the approach in terms of 
chaotic transients, like in (Lagha & Manneville, 2007a). Its detection via a single experiment 
where R was progressively decreased by small steps has been judged sufficient for the present 
purpose. 

Finally, in a very wide domain 680x340 of size comparable to that of the largest experimental 
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Figure 4: Patterning in a domain of size 128x84. Left: Patterns at R = 154. Right: Bifurcation 
diagram (distance to laminar flow as a function of Reynolds number, see text). 



Figure 5: Snapshot of the solution in a domain of size 680 x 340 at R — 151. To reach such 
a size, the in-plane resolution has been lowered to N x = L x , N z = L Z) without destroying the 
pattern. 


setups (Prigent et al. 2002), patterns with many wavelengths were obtained. Comparing £gure[5] 
and figure 8 in (Manneville & Rolland 2011) one note that the model generates outputs quite 
similar to what is obtained in DNSs at reduced wall-normal resolution, itself representing the 
experimental situation reasonably well. 

Whereas at a qualitative level the model is fully satisfactory, we have however noted certain 
quantitative discrepancies. First the transitional range is shifted downwards somewhat more 
importantly than in the DNSs (Manneville & Rolland 2011). This can be understood by noticing 
that the amount of energy extracted from the base flow by viscous stresses at the plates is 
transferred though a very short range of wall-normal small scales most likely to dissipate it. This 
artificially maintains more turbulent activity in the system at given R , or equivalently a similar 
turbulence level at lower R , than in better resolved simulations and laboratory experiments 
where energy is transferred to and efficiently dissipated in much smaller scales. (There is 
little or no trade-off for the in-plane dissipation that is treated like in the full-3D DNSs.) 
For a concrete comparison, experiments (Prigent et al. 2002) and well-resolved full-3D DNSs 
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(Duguet et al. 2010) give the upper threshold (featureless turbulence) at R t ~ 410 and the 
lower threshold (global stability) at R g ~ 325. In simulations at reduced resolution (Manneville 
& Rolland 2011), N y being the number of Chebyshev polynomials used in the representation 
of the wall-normal dependence, these values are shifted down to 350 and 270 for N y = 15 and 
to 275 and 215 for N y = 11. Here we have R t ~ 159 and R g ~ 150 but the pattern is still 
well rendered. This larger shift can therefore be understood because the effective wall-normal 
resolution is much lower. The fact that a physically relevant solution is obtained here while the 
Chebyshev implementation breaks down with similarly few modes is due the optimal rendering 
of boundary conditions on v achieved by our basis choice (see below). 

Second, though the angle between the bands and the streamwise direction is correct, the 
wavelengths of the pattern, both streamwise and spanwise, are too short by a factor of 1.5 to 2, 
and the pattern’s orientation in domains of intermediate size fluctuates more than in the DNSs. 
The amount of enhancement is however difficult to appreciate quantitatively. These phenomena 
remain unexplained for the moment but might relate to the effect of the wall-normal resolution 
on the streamwise coherence that was shown to play an important role on the existence of the 
pattern (Philip & Manneville 2011). A hand-waving confirmation of this effect on the robustness 
of the bands comes from the continuous trend observed as the resolution is decreased, here as 
the truncation level i max is lowered, in rough correspondence with what was observed when 
reducing the wall-normal resolution in DNSs. First, a conspicuous steady pattern is observed 
with i m ax = 3. Next, for i max = 2 (not reported here but studied in parallel) coexistence of 
fluctuating, wide, laminar and turbulent domains are observed in an even narrower Reynolds 
range; these domains remain disorganised and do not form bands. Finally, for i max = 1 (Lagha 
& Manneville 2007a), streaks stay short, the transitional range seems to be reduced to a point 
at a somewhat larger value ( R g ~ 170), and wide steady coexisting domains do not exist: either 
turbulent domains grow from small germs in a laminar background or the reverse (Manneville 
2009). As a matter of fact, a similar but worse situation occurred in full 3D simulations at 
exaggeratedly reduced resolution since nontrivial states with unphysical small scales and no 
patterns were obtained for N y = 9 and 7 and blow-up occurred for N y < 7, which give a 
marked advantage to our separately optimal wall-normal representations of v and (. These 
observations should contain some physics that warrant to be elucidated, as suggested in the 
next section. 


4 Perspectives 

Understanding the transition to turbulence in wall-bounded flows, and especially PCF that is 
linearly stable for all R and displays alternating laminar and turbulent oblique stripes on its 
way to fully developed turbulence, is a hard problem when starting from the NSE. Some sim¬ 
plification can be expected by taking a key ingredient into account: the transition takes place 
at moderate values of the Reynolds number for which the flow is controlled by the presence 
of coherent structures (Hof al al. 2004; Bottin et al. 1998). The model that we have derived 
incorporates this feature by means of a Galerkin expansion of the dynamics optimally adapted 
to the boundary conditions at the walls. Truncating the expansion beyond lowest nontrivial 
order, keeping 7 amplitudes instead of 3 in (Lagha & Manneville 2007a,b), has allowed us to 
recover the experimentally and numerically observed patterning at minimal price (downward 
shift of the transitional range, somewhat too short pattern wavelengths). This limitation can 
be overcome by increasing ?' max , which raises the interesting question of the rate of conver- 
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gence of the approximation. Such a study would possibly be rewarding because, to be precise, 
DNSs treating the three space directions on a similar footing are computationally extremely 
demanding (Duguet et al. 2010). In a spirit akin to that of large eddy simulations, the present 
modelling avoids to waste numerical resources by singling out hydrodynamic coherence in the 
wall-normal direction and skipping the explicit computation of the small scales. 

An indication that the convergence of our approach could be fast is the observation reported 
by one of us - see Fig. B6 in (Mannevillc 2015) - that the fraction of the perturbation energy 
retained in projections of fully resolved numerical solutions onto the basis considered here 
tends to 1 exponentially fast as the truncation order is increased: 90%, 97%, and 99% for 
imax = 1, 3, and 5, respectively. This does not prove that the global dynamics of the system 
would be equally well captured quantitatively by increasing the number 2f max + 1 of fields 
in the model but hints at such a convergence, as generally expected for spectral approaches, 
here relying on specific complete series of Jacobi polynomials (Rolland 2012). At the price of 
a pre-treatment of the problem that amounts to the once-for-all automated derivation of an 
effective set of equations of sufficiently high order, it might be found interesting to replace the 
full 3 D numerical simulations of the NSE by a finite set of two-dimensional partial differential 
equations already taking the continuity condition fully into account and managing with wall- 
normal coherence in the transitional range. A quantitative estimate of the expected gain in 
terms of memory requirements and time steps definitely warrants further study. 

In a complementary perspective, one can rather think of analysing the properties of the 
model. First, in-plane coherence may be added to the wall-normal coherence inherent in the 
derivation. This can be done by inserting specific assumptions about the (x, ^-dependence of 
fields 0 and 0, in particular strict periodicity in space at the scale of the MFU (Jimenez & 
Moin 1991). With z max = 1 and further limiting the in-plane expansion to the first harmonic, it 
is then straightforward to recover Waleffe’s models (Waleffe 1997) by making the corresponding 
educated guess. A system of eight equations for eight amplitudes is obtained, identical to his 
system (10) but with a different set of coefficients acknowledging the difference in boundary 
conditions (which, in passing, shows the structural genericity of that model). For example, the 
equation for the streamwise mean-flow component called M by Waleffe and governed by his 
equation (10a) here reads: 

ipi - ypuUx = \rsi 01 [(a 2 + 7 2 )BE - 2 UV] 

where U\ = M— 1, while other symbols have the same definition as in (Waleffe 1997), especially 
the streamwise and spanwise wave-vectors a = 2 tt / i x and 7 = 2ttH z , i x and l z being the 
dimensions of the MFU. The numerical values of coefficients in the equation above can be 
obtained from the formal expressions in [A] Following the very same line, a study to be presented 
elsewhere (Mannevillc, in preparation) shows that uniform large scale flows are generated just 
by shifting the phase of specific ingredients of Waleffe’s eight-equation model. Combining this 
to the introduction of appropriately weighted in-plane second harmonics should help us to 
account for oblique coherent structures like those recently found by Daly and Schneider (2014), 
though the actual derivation of a model possessing them as fixed points would be cumbersome. 

Beyond the simple hypotheses corresponding to strictly periodic coherent structures, the 
next step is to describe spatially slow turbulence modulations corresponding to the patterns 
observed experimentally through the formal introduction of a slow dependence of the amplitude 
of the local bifurcated state, in the spirit of the derivation of standard multiple-scale envelope 
equations. The approach cannot be made as rigorous as, e.g., for convection the since the 
bifurcated state stays at finite distance from the laminar-flow base, which leaves room for 
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further modelling. Of the two scales introduced, the fast one accounts for mechanisms at the 
MFU scale and the slow one corresponds to the modulations. The slow variables are driven 
by source terms arising from a filtering of the Reynolds stresses, like in (Lagha & Manneville 
2007b) and it can be seen that the modulation of the uniform large scale flows alluded to 
above generates nonlocal contributions of the class identified by Hayot and Porneau (1994) as 
playing an essential role in the balance between laminar and turbulent regions responsible for 
patterning. But, in contrast with their phenomenological introduction of such contributions, 
here they directly arise from the equations and are therefore sensitive to the local orientation 
of the flow with respect to the streamwise direction, hopefully giving a microscopic support to 
the empirical observations of Duguet and Schlatter (2013). 

Finally, large scale flows are present already with i max = 1 (Lagha & Manneville 2007b; 
Manneville, in preparation), though steady patterns are not observed in that case (Manneville 
2009). Taking smaller wall-normal scales into account (i max > 1) is therefore necessary for a 
theoretical interpretation of the stabilisation of long-wave modulations observed with i max = 3, 
as reported in (j3j Simplification of models with higher truncation levels would then take 
advantage of the diagonal dominance of matrices A, P, and P noticed earlier to perform the 
adiabatic elimination of terms of least relevance yielding an effective model for the slowly 
evolving terms. Such a heavy work could however possibly not be necessary and considering 
seven fields might be sufficient up to an optimisation of the model’s coefficients. As a matter of 
fact, the three first amplitudes (To, Ti, dq) are the most appropriate to deal with the nontrivial 
properties of the in-plane flow dependence. So, if one is willing to include more of the wall- 
normal dependence, it should suffice to consider that the pairs (T 2 , T 2 ) and (T 3 ,<P 3 ) collect 
all the higher order contributions of each parity and, owing to its generic structure, to restrict 
oneself to the consideration of the seven-held model as an effective system replacing the NSE. In 
this perspective, except as a starting guess, sticking to the values of the coefficients obtained in 
the strict Galerkin expansion is not advisable and introducing some multiplicative randomness 
at appropriate strategic places like in (Barkley 2011b) seems profitable. Applying the program 
sketched in the previous paragraph to this new primitive problem is currently developed, which 
is expected to improve over the one-dimensional phenomenological approaches of Manneville 
(2012) and Hayot & Porneau (1994). 

5 Conclusion 

The subcritical coexistence of different regimes forming laminar and turbulent patterns in PCF 
and other wall-bounded how configurations is a difficult problem in which the interplay of 
mean how corrections and finite amplitude perturbations plays a crucial role. Our approach 
via Galerkin decomposition yields explicit models replacing the NSE by coupled systems gov¬ 
erning amplitudes that encode the gross features of the how. The derivation is systematic 
and the structure of the obtained models is generic, reflecting that of the primitive equations. 
Simulations of those models reproduce the patterning provided that the truncation level is 
not too low. They are next amenable to further analysis, especially through in-plane space 
dependence assumptions and explicit scale separation. Here, this program has been developed 
for PCF but its adaptation to other hows such as plane Poiseuille or Couette-Poiscuille how, 
cylindrical Couette-Taylor how, etc. is straightforward. Obtaining a Barkley-like model for 
Hagen-Poiseuille how from hrst principles can also be considered along similar lines, using ba¬ 
sis functions adapted to the tube geometry and the no-slip condition. The extension to the 
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less trivial case of boundary layer flows of various kinds, with their free-stream boundaries 
at infinity, remains a stimulating challenge. Once obtained such models offer tools to scruti¬ 
nise laminar-turbulent coexistence and provide us with detailed physical explanations of this 
phenomenon of great conceptual and practical importance. 


A Explicit expressions 

A.l Basis functions 

• In-plane velocity components and wall-normal vorticity component: 

/(0) = |A5(1 -y 2 ), 

/(l) = JVT05 (1 - 2/ 2 ) 2/, 

m = f v'h (l - </ 2 ) (</ 2 - |) , 

/(3) = f \/Tl55 (l — y 2 ) (y 2 — §)?/, 

/(4) = §Vfm(l-y 2 )(y 4 -^y 2 + ±), 

f (5) = fv / 70(l- 2 / 2 )(i/ 4 -i 2 / 2 + 1 i)|/,... 

• Wall-normal velocity component: 

9(1) = r e VS5(l -y 2 ) 2 , 

g( 2) = ^a/385 (l-y 2 ) 2 y, 

g ( 3 ) = 1>/91 (t — 2 / 2 ) 2 ( 2 / 2 — n ) , 

<7(4) = 1V385 (y 2 -l) 2 ( y 2 -^)y, 

9( 5) = Ill's/1309 (y 2 — l) 2 (y A — \y 2 + ^) , ■ ■ ■ 


A.2 Coefficients in the evolution equations 


Taking care of the order of the subscripts introduced, symmetries within the sets of coefficients 
are easily detected, directly or via integration by parts. Energy conservation relies on the 
symmetries of coefficients introduced in the expressions of nonlinear terms. The elements of 


matrix C appearing in (8) are straightforwardly obtained as Cji = dy fjg[. 


A. 2.1 Equation (J 5 j) for 


• linear terms: 

matrix I: 8# = dy g 3 (; yg t ), 

matrices A and A: a 3% = f^' dyg 3 g”, a 3i = /+/ dyg 3 (yg"), 
matrix P: p jt = dy g 3 g’l"\ 

• nonlinear terms, for j e (1 : i max ): 

Imax 7.max ^max 7max 

= EE Qjik A [d x (U i Vk) + d z (W i V k )] + EE Qjik^(l^Vk) 

i =0 k= 1 i= 1 k= 1 
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-max ‘max 


r i ik (^) + 2d ™ + d ** ( W i W *)\ 

i= 0 fc=0 

*max *max 

-EE f jik [d x {U t V k ) + d z {Wi\ 4)], (10) 

i=0 k=1 

with: q jik = /J^ 1 d yg 3 fig k , qpk = C! dyg^g^y, 

r jik = Cl dy gj ( fifk)', f jik = Cl d y <g 3 ( flg k )". 

A.2.2 Equation (|g]) for Tj 


• linear terms: matrices B, B, and P: 

bji = C! d y fj ( yfi ), hi = Cl dy fai, Pji = Cl dy flf", 

• nonlinear terms, for j G (0 : ?' max ): 

^max *max 

n} z) = w i w k) + ( d zz- d xx)(Ui w k )\ 

i= 0 k= 0 

^max *max 

+ EE s jik [d z {UiV k ) - d x {W.y k ))l (11) 

i= 0 fc=l 

with s jik = Ci dy fjfifk, Sjik = Cl d y fjifidk)'- 
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